#ifndef THEOR_BI_H
#define THEOR_BI_H

#include "commonheader.h"
#include "gslheader.h"
#include "int_1d.h"
#include "spline_bessel.h"

#define Z_SPACING 0.25
#define XMIN 12500.0
#define XMID1 13700.0
#define XMID2 14300.0
#define XMAX 15500.0
#define PTS_MIN 96
#define PTS_MID 96
#define PTS_MAX 96
#define NSUB_MIN 1
#define NSUB_MID 1
#define NSUB_MAX 1
#define PTS_KINT 96
#define NSUB_KINT 200
#define PTS_CL 96
#define NSUB_CL 500
// These constants should, in a future version, be determined dynamically
// from the transfer function.

typedef struct{
  int iso[2], pol[2];
  int li[3];
  double x, qmin, qmax, power;
} funcparams;

typedef struct{
  int num_shape, num_iso, num_pol;
  int num_shape_do, num_iso_do, num_pol_do;
  int do_iso_summ, do_pol_summ;
  int *shape_arr, **iso_arr, **pol_arr;
  int *iso_summ, *iso_summ_nl, *pol_summ, *pol_summ_nl;
  double *tf_corr, **red_tilt, **rel_ampl;
  double pivot, red_tilt_adi, red_tilt_iso, rel_ampl_ai, rel_ampl_iso;
  int num_l, num_p;
  int *l_arr;
  double *p_arr;
  gsl_spline **spline_arr;
  gsl_interp_accel **acc_arr;
} transferfunc_vars;

typedef struct{
  int num_li;
  int num_x;
  int *li_arr;
  double *x_arr;
  double *w_arr;
  double *alpha_arr;
  double *beta_adi_arr;
  double *beta_iso_arr;
  double *gamma_arr;
  double *delta_arr;
} k_integral_vars;

double *read_transfer(char filename[], int first);
void make_splines(double *ptf_arr);
void free_splines(void);
void init_bess_table(double x_max, double z_spacing);
void init_theor_bispec(char **filename);
void free_all_splines(void);
double ki_integrand(double k, void *pparams);
double tf_integrand(double k, void *pparams);
void init_kintvars_x(void);
void init_kintvars_w(void);
int adaptive_step(int li, int stepstep, int stepmax);
void init_kintvars_l(int limin, int limax, int stepstep, int stepmax,
		     int do_shapes);
void calculate_k_integrals(funcparams *pparams, int do_shapes);
double *compute_bispectrum_local(funcparams *pparams, int shape);
double *compute_bispectrum_local_indep(funcparams *pparams, int shape);
double *compute_bispectrum_equil(funcparams *pparams, int shape);
double *compute_bispectrum_ortho(funcparams *pparams, int shape);
double *compute_cl(funcparams *pparams);

#ifndef THEOR_BI_GLOBAL_VARIABLES

extern transferfunc_vars tfvars;
extern k_integral_vars kintvars;

#endif

#endif
